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AUTOMATED PROCEDURE FOR DESIGN OF WING STRUCTURES 
TO SATISFY STRENGTH AND FLUTTER REQUIREMENTS 

By Raphael T. Haftka* 

Langley Research Center 

SUMMARY ‘ ' 

A pilot computer program has been developed for the design of minimum mass wing 
structures under flutter, strength, and minimum gage constraints. The wing structure is 
idealized by finite elements, and second- order piston theory aerodynamics is used in the - 
flutter calculation. Mathematical programing methods are used for the optimization. 

Computation times during the design process are reduced by three techniques. 

First, iterative analysis methods are used to reduce significantly reanalysis times com- 
pared with the original analysis of the structure. Second, the number of design variables 
is kept small by not using a one-to-one correspondence between finite elements and design 
variables. Third, a technique for using approximate second derivatives with Newton's 
method for the optimization is incorporated since it proves to be superior to the commonly 
used Davidon-Fletcher-Powell algorithm. 

The program output is compared with previous published results. Reasonable agree- 
ment is obtained. In addition, it is found that some flutter characteristics, such as the 
flutter speed, can display discontinuous dependence on the design variables (which are the 
thicknesses of the structural elements). It is concluded that it is undesirable to use such 
quantities in the formulation of the flutter constraint. 

INTRODUCTION 

The usual design process for aircraft wing structures consists of sizing, and per- 
haps optimizing, for strength (with the use, for example, of fully stressed design tech- 
niques, ref. 1), checking for flutter and, if required, determining a flutter fix by trial and 
error. In the past few years, there has been considerable interest in developing optimi- 
zation procedures for the design of aircraft structures to satisfy aeroelastic constraints. 

Work in this field is reviewed in references 2 and 3. When flutter is the only con- 
straint (except for minimum gage) optimality criteria (refs. 4 to 6) or special flutter 
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oriented optimization algorithms (ref. 7) can be used. When flutter is only one of many 
design conditions, the more general mathematical programing methods are more con- 
venient. The first effort in this area, dealing with the design of an airfoil for minimum 
wing drag work during a given flight mission, is presented in reference 8 where flutter 
Mach number is one of several constraints included. A gradient method is used to opti- 
mize the two design variables which are the thickness and chord of the airfoil. 

An automated design program called SWIFT is presented in reference 9, wherein a 
simple structural model of an aircraft wing is used to show the effects of strength and 
flutter requirements on the design of minimum mass aircraft wing structures. The wing 
is idealized as an isotropic sandwich plate with a specified variable depth between covers 
and a variable cover thickness distribution which is optimized. The sequence of uncon- 
strained minimizations technique (SUMT, refs. 10 and 11) with an interior penalty function 
is used. The Davidon- Fletcher -Powell algorithm (ref. 12) is used for each unconstrained 
minimization. Reference 13 is similar to reference 9 in that it also uses a plate wing 
model but differs from reference 9 in that it employs the feasible directions method 
(ref. 14) to optimize both the structure and the configuration under a large number of 
constraints. 

The plate model is not adequate for general built-up wings, for which finite element 
models are required; such models were used even in some of the earlier works (refs. 3, 

4, 6, and 7). Finite element analysis, however, involves considerably more computation 
than plate analysis. In the past when finite element models have been used in conjunction 
with mathematical programing methods, computation times tended to become excessive. 
For example, the optimization of a finite element wing model under strength, buckling, 
frequency, and flutter constraints is discussed in reference 15 with the use of the same 
optimization algorithm as in reference 9. Computation times are long even for crude 
finite element models. 

The present paper considers the structural optimization for minimum mass of built- 
up wings (ribs, spars, cover panels) under maximum stress, minimum gage, and flutter 
constraints. Particular attention is given to the problem of reducing the computing time 
required for executing hundreds of analyses of a finite element wing model, as required by 
mathematical programing techniques. Of course, computer time is not the only consider- 
ation in computer cost, but the other considerations, such as core storage, are installation 
dependent and therefore not discussed herein. It is shown that the use of appropriate 
analysis algorithms can reduce the computation time of a reanalysis of the wing to a small 
fraction of the time needed for the original analysis. In addition, a technique for using 
approximate second derivatives in conjunction with Newton's method is shown to be 
superior to the Davidon-Fletcher-Powell (DFP) algorithm. Attention is also given to the 
efficient choice of design variables that will reduce their total number. 
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A computer program using, these ideas was written for the design of wings assumed 
to have symmetric cross-section profiles and to be clamped at the root. The flutter cal- 
culations are made at a fixed altitude with the use of second-order piston theory aerody- 
namics and natural vibration modes as generalized coordinates. Both Newton’s method 
and the DFP algorithm were used in the present program which is called WIDOW AC 
(Wing Design Optimization With Aeroelastic Constraints). Analysis and sample results 
are presented along with a study of discontinuous behavior of the flutter speed as a func- 
tion of design variables, which was found to occur under some circumstances. 

SYMBOLS 

4 

The physical quantities used in this paper are given in both the International System 
of Units (SI) (ref. 16) and in the U.S. Customary System of Units. The measurements and 
calculations were made in U.S. Customary Units. Appendix A presents some factors 
relating these two systems of units. 

Ajj generalized aerodynamic forces (see eq. (13)) 

a speed of sound 

a l’ a 2’ a 3 modal amplitudes 

[B] matrix defined by equation (21) 

f(v) object function for optimization 

G shear modulus 

g structural damping 

gj structural damping associated with jth vibration mode 

g l’ g 2i’ g 3i’ g 4j constraint functions 
h total local depth of wing 

[I] unit matrix 

[K] stiffness matrix 
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load vector 


(L) 

M Mach number 

[M] modal mass matrix 

n number of vibration modes 

P(r) sum of object function and penalty terms 

Ap * amplitude of lifting pressure (see eq. (15)) 

Apj Ap due to oscillation of wing in'jth vibration mode 

Ap(x,y,t) lifting pressure at point x,y on wing at time t, positive up 

r factor of penalty function (see eq. (5)) 

S area 

[S] stress matrix 

(s) stress vector 

t time 

t e j prescribed minimum thickness for jth element 

ti thickness at ith node 

tj thickness of jth element 

ts,min minimum skin thickness 

V free-stream velocity 

Vf flutter speed 

Vf cr critical flutter speed 
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v vector of design variables 

Vj^ ith design variable 

W amplitude of normal displacement (see eq. (15)) 

W normal displacement at wing surface 

W for oscillation of wing in ith vibration mode 
(w) displacement vector 

x,y,z chordwise, spanwise, and normal coordinates, respectively 

y ratio of specific heats 

Kronecker delta 
v Poisson’s ratio 

p air mass density 

i 

a cr critical value of stress 

o x ,Oy,T X y components of normal and shearing stress in wing elements 

co frequency, rad/s 

coj flutter frequency 

coj frequency of jth vibration mode 

ANALYTICAL PROCEDURES 
Problem Definition 

The general optimization procedure is to determine a vector v that solves the 
problem: 

Minimize the object function f(v) subject to the constants gj.(v) =0 (j = l(l)m) (1) 
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Herein the object function is the structural mass of the wing, and the following three con- 
straints are considered: 

(1) Flutter constraint: The flutter speed Vf at a specified altitude is greater than 
a prescribed value Vf cr . Thus, 


ei s 1 ' 



( 2 ) 


(2). Stress constraints: The Von Mises stress a at any point in the structure does 
not exceed a prescribed value o cr . Thus, 


&2i = 1 


2 0 


cr,i 


(i = l(l)n el ) (3) 


where and a cr j are the Von Mises stress and the critical stress for the ith ele- 
ment, respectively, and n e ^ is the number of finite elements in the structure. 

(3) Minimum gage constraints: No element thickness tj shall be less than a pre- 
scribed value t e j. Also, at no node on the cover shall the skin thickness t^ be less 
than a prescribed value t S inin . (The double constraint is needed because, on the cover 
element, thickness is the average of nodal thicknesses and a constraint on the element is 
not sufficient.) Thus, 


&3i 


= l ts > min > Q 
*1 


,. 1 - 41*0 


g4 i y 


( l = ^Kode) 

(j = l(l)n el ) 


( 4 ) 


Optimization Technique 

■ The mass is minimized subject to these three constraints by using the SUMT 
approach (refs. 10 and 11) with an interior penalty function. With this method the problem 
defined by equation (1) is replaced by 

Find lim min (P(r)) where 
r— 0 

m 

P( r ) = f(v) + r J -±- (5) 

i=l ^ v > 


For each value of r an unconstrained minimization of P is carried out. 
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The DFP algorithm is among the most efficient methods for unconstrained optimi- 
zation which are commonly used with the SUMT approach (refs. 9 and 15). Such methods 
require calculation of the first derivatives of the object function, and each unconstrained 
minimization is executed by about n one-dimensional searches, n being the number 
of design variables. For this investigation, however, Newton's method, which requires 
both first and second derivatives of the object function, is found to be more efficient than 
the DFP algorithms. It is shown in the following paragraphs that the first derivatives of 
the constraints (which are needed anyhow for calculation of the first derivatives of the 
P function) are sufficient to establish a good approximation to the second derivatives of 
the P function. With Newton's method, each unconstrained optimization requires only a 
small number of one -dimensional searches, a number which is independent of n.' The 
second derivatives of P are obtained as follows: 

Differentiating P twice, one has 


a 2 f 


9 2 P = 

8Vj 9Vj 9Vj 9Vj 


V L ag i 9g i a2g i \ 

+ r Zfe^" g ia^j g i 


( 6 ) 


An approximation to 


9 2 P 
9vj 9Vj 


is available under the following assumptions: 


(1) known 

(2) Some of the constraints are nearly critical so that for some values of 

i, g^v) is very small 

(3) The second derivatives of the g's are not very large 

Under these assumptions, the second term in the summation in equation (6) may be 
neglected so that 


9 2 P s 8 2 f ( 0 

9Vj 9v, 9V; 9vj jl 


(7) 


,j „ V J 


where 


11 4 \ 9v j 8v ij 1 


In the present work, 


9 2 f 


9Vj 9Vj 


= 0 because the mass is a linear function of the design 


variables, and the first derivatives of the g's are obtained numerically by forward 
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differences (with reanalysis of the structure). The matrix Q is positive semidefinite, 
but when the number of active constraints is small, it is singular or ill-conditioned. This 
is because when there are only a few constraints and they are linear with respect to the 
design variables, there is no unique optimum point (the nonlinearity of the constraints 
may create a unique optimum point). To guard against this possibility, Qjj is replaced 

by 


Q ij - Q ij C 1 + e6 ij) 

A value -of the constant e of 0.01 has been found satisfactory to avoid the singularity or 
ill - conditioning. 


Idealization of Structure and Design Variables 

The wing is assumed to have a symmetric airfoil profile and to be clamped at the 
root. Constant strain triangular membrane elements are used to represent the wing 
cover panel, and quadrilateral shear-web elements represent the ribs and spars (see 
appendix B for additional details). 

The natural choice for design variables in a structural design problem is the 
dimensions of the structural elements. When the analysis of the structure is carried out 
by using a finite element idealization, the thicknesses of the finite elements are often 
picked as design variables. This choice, however, is not always logical. For example, 
an area of stress concentration may require a fine mesh of finite elements, but it is not 
always reasonable to assign to the thickness of that area a large number of design vari- 
ables. When the optimization is done by mathematical programing, it is often necessary 
to use far fewer design variables than finite elements and, hence, distinguish the finite- 
element model from the design variable model (ref. 17, for example). Typical finite ele- 
ment and design models are shown in figure 1. In WIDOW AC, the design variables do not 
describe the thickness of individual finite elements but of segments of the wing, each seg- 
ment containing several finite elements. The wing planform is divided into a number of 
triangular or quadrilateral segments. The division is such that no finite element crosses 
the boundary of a segment. The design variables are the thicknesses of the cover panels 
at the vertices of the segments (indicated by the numbers in fig. 1). The thickness varies 
linearly in the triangular segments and is piecewise linear inside the quadrilateral ones 
(each quadrilateral is divided into two triangles and the thickness variation in each one is 
linear). It is optional to lump together several vertices and assign to them the same 
design variable. For example, if the same design variable is assigned to vertices 30, 37, 
38, and 33 in figure 1, segment VI will be a constant-thickness segment. Similarly, all 
the shaded segments in the figure can be lumped together. On the other hand, it is 
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optional to create discontinuities in thickness by assigning different design variables to 
the same point, if that point is a vertex of more than one segment. 

The shear webs which represent the spars and ribs are assigned design variables by 
finite elements, but it is again possible to lump together several elements. For example, 
the thicknesses of all the spar elements may be lumped together as one design variable 
and the thicknesses of all the rib elements as another. 

The methods just discussed for assigning design variables provide considerable 
flexibility. It is possible, for example, to first optimize the wing with a low number of 
design variables. This first optimization usually discloses that certain regions of the 
wing are more critical than others. It is possible as a next step to increase the number 
of design variables that describe these regions but continue with a crude representation of 
less important regions. 

Two basic questions have to be answered regarding the usefulness of this method 
for reducing the number of the design variables. First, what is the difference in the final 
mass because of the smaller number of design variables? Second, how different are the 
resulting designs ? To be useful for preliminary design it is important that basically the 
same design will emerge with the smaller number of design variables as with the larger 
one. These two questions have been investigated for the example problems and the 
results, presented in a subsequent section, are encouraging. 

Flutter and Stress Analyses 

The calculation of the flutter speed and stresses appearing in constraint equa- 
tions (2) and (3) is given as follows. The displacement method of finite element analysis 
is used so that the displacements are the primary variables. The displacement vec- 
tor . (w) is found from the following equation: 

[K](w} = (L) (8) 

where [K] is the stiffness matrix and {l} is a load vector. Stresses are calculated 
from the displacements for each element 

(s) = [S](w} (9) 

where [S] is the stress matrix and (s) is the stress vector which has the components 
ct x , <7y, and T xy for each element. The Von Mises stress a which is used in con- 
straint equation (3) is 
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(10) 



°x°y 


+ 3 T xy 2 ) 


1/2 


For the flutter problem the load vector {L} in equation (8) is replaced by inertia, 
structural damping, and aerodynamic forces. When harmonic motion is assumed, the 
flutter equation is 


[[K] (1 + i[G]) - u> 2 [M] - [A(a),V)]] (w) = 0 (11) 

where [A] is the aerodynamic matrix, which (for a fixed altitude) is a function of the 
frequency o> and speed V, [G] is the structural damping matrix, and [M] is the 
mass matrix. Equation (11) is an eigenvalue problem in the pair V,u> determined from 
the condition that the determinant of equation (11) should vanish. 

To reduce the order of the problem it is common to use the natural vibration modes 
as generalized coordinates for the flutter problem. In terms of the natural vibration 
modes (normalized to give a unit generalized mass matrix) the determinant of equa- 
tion (11) is 


A »i 


(oo* 






: j) 


6 ij 


= 0 


(i,j = l(l)n) (12) 


where coj is the jth natural frequency, gj is the structural damping associated with the 
jth mode, and A^j is the generalized aerodynamic force 

A ij = ^Ap^dS (13) 

S 


The pressure loading on the wing as predicted by second-order piston theory (refs. 18 
and 19) is 


Ap(x,y,t) = -2pa 1 + 


^ + -± M + v ‘ W 

9x/ \9t 9 x) 


(14) 


When harmonic dependence on time is assumed so that 
Ap = Ape 

> 

W = We iwt 


(15) 


equation (14) becomes 
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( 16 ) 


Ap = -2pa jl + M — ViwW + V —■ 


Equation (16) is substituted into equation (13) to calculate the generalized aerodynamic 
forces. 


Solution Methods 


The choice of computation algorithms for a design program is not motivated so 
much by their efficiency in the execution of one analysis but by their efficiency in exe- 
cuting many analyses. Specifically, iterative methods which may not be appropriate for 
a single analysis because they must be started with a very good guess of the solution 
become very attractive when many repeated analyses are to be executed. The choice of 
appropriate computation algorithms, described in the following sections, in WIDOW AC has 
reduced the computer time required for a reanalysis to a small fraction of the computer 
time required for the first analysis. The following steps are required for the computa- 
tion of stresses and flutter speed: 

(1) Stiffness matrix assembly 

(2) Solution for displacement and stresses (eqs. (8) to (10)) 

(3) Calculation of vibration modes 

(4) Flutter solution 

The algorithms used in these steps are discussed in the following sections. 

Stiffness matrix assembly .- The design variables describe the thicknesses of the 
structural elements, and because no bending elements are used, the stiffness matrix 
depends linearly on the design variables. This linear dependence is used to reduce the 
assembly time during the optimization process by computing the stiffness matrix as 


[K] = [K] q + 


9K 

3V; 


i=l 


(17) 


where [K] is the stiffness matrix, v^ is the ith design variable, and [K] 0 is the 
residual stiffness matrix obtained when all the design variables are set to zero (i.e., the 
stiffness of the part of the structure which is not optimized). It should be noted that, 
although higher order finite elements are not presently available in WIDOWAC, this 
method of assembly permits the use of such elements with almost no penalty in computa- 
tion time because the element matrices are computed only once. 
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Solution for displacements and stresses .- The matrix equations resulting from the 
finite element model are sparse. Efficient handling of operations with sparse matrices 
is achieved by using the node as the basic unit of the calculation rather than the degree of 
freedom. Each node has three degrees of freedom (translations in the x, y, and z 
directions), and this situation results in a stiffness matrix built up of 3 x 3 submatrices. 

The use of 3 x 3 submatrices has two advantages. First, the storage requirements 
and amount of bookkeeping operations required for taking full advantage of the topology of 
a sparse matrix are greatly reduced. Second, operations on 3 x 3 matrices are programed 
without use of loops by writing out all the steps in the operations. As a result computa- 
tion times with these matrices are reduced by a factor of 3 or 4. 

Calculation of vibration modes .- The first step in the flutter analysis is the calcula- 
tion of a number of the lowest vibration frequencies and the corresponding vibration 
modes. The eigenvalue problem has the following form: 


[[K] - u> 2 [M]] (w) = (0) (18) 

By introducing 

(w) = [M] 1/2 (w} (19) 

the eigenvalue problem (eq. (18)) may be transformed into 



1/2 


[K] _ 1 [M] 


.1/2 _ J_ 


0J* 


[I] 


(w}=(0) 


( 20 ) 


so that the problem is reduced to finding the largest eigenvalues and corresponding eigen- 
vectors of the matrix [B] 

[B] = [M] 1 / 2 [K]' 1 [M] 1/2 (21) 

The eigenvalue problem is solved by a simultaneous iteration method (ref. 20) coded 
by Alfred F. Vachris, Jr., of the Grumman Aerospace Corporation. The essential feature 
of the method is that it is very efficient when a good approximation to the eigenvectors is 
available. If the relative error in the eigenvalue is O(e) with e < 1, it is reduced in one 
iteration to 0(e2). The method is, therefore, very suitable for the optimization process. 

Solution of flutter determinant .- The flutter determinant (eq. (12)) is a function of V 
and cu. The usual way of finding values of V and u> which make the determinant van- 
ish is to evaluate all eigenvalues as a function of the reduced frequency (ref. 21) leading to 
the familiar V-g diagram method. 
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In WIDOW AC, a different method is employed which is more efficient when a very 
good initial approximation to the solution is available. The method (same as ref. 9) uses 
a Newton iteration with V and o> to drive the real and imaginary parts of determinant 
to zero simultaneously. For simplicity, the calculation is done with zero structural 
damping gj and for a constant air density (i.e., constant altitude). 

For each step in the course of the optimization, the initial approximation for V 
and o» is based on extrapolation of the results from preceding steps. Only one of the 
flutter modes is followed and therefore there is a danger that during the design process 
a different mode will have a lower flutter speed. Because of this possibility, the final 
design should be checked for a wide range of V and u>. 

The flutter point for the initial and final design are calculated by separate runs, in 
which the flutter determinant is computed at a grid of points in the V-u> plane, and 
points are found where both the real and imaginary parts of the determinant change sign. 
The method may be very inconvenient and inefficient for one analysis, but it is very 
efficient for design. 


Computer Program 

A pilot computer program has been developed for the design of minimum mass wing 
structures under flutter, strength, and minimum gage constraints. The wings are 
assumed to have symmetric cross-section profiles and to be clamped at the root. The 
flutter calculations are made at a fixed altitude with the use of second-order piston 
theory aerodynamics and natural vibration modes as generalized coordinates. Both 
Newton’s method and the DFP algorithm are used in the present program which is called 
WIDOW AC (Wing Design Optimization With Aeroelastic Constraints). The pilot computer 
program was developed for use on the Control Data Series 6000 computer systems. 

APPLICATIONS AND RESULTS 
Full- Depth Sandwich Wing 

Results obtained from WIDOW AC for the titanium clipped delta wing, defined in 
figure 2, are compared with similar results obtained from the SWIFT program of ref- 
erence 9. SWIFT is a wing design program which served as the starting point for 
WIDOW AC development. It treats a wing as a sandwich plate with a core having infinite 
transverse shear stiffness. Plate theory is used for the structural representation. The 
design variables are the coefficients of a polynomial which expresses the thickness of the 
panel covers over the wing. In the numerical analysis, the spanwise behavior is dis- 
cretized by a finite difference scheme; the chordwise behavior is handled by a parabolic 
approximation. 
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.The SWIFT model had seven finite difference stations, and this modeling results 
effectively in 21 degrees of freedom. The finite element model used in WIDOW AC must 
consist of many more degrees of freedom than a plate wing model used in SWIFT to pro- 
vide adequate accuracy; a model with 93 degrees of freedom was used in the calculations 
for the present paper. To simulate the rigid core of the SWIFT model, the shear webs 
in WIDOW AC were made very stiff and were not assigned any design variables. Only the 
thickness of the cover panels was optimized. 

The division of the wing into segments is shown in figure 1. The optimization 
process was usually started by assigning one design variable to each segment. In this 
application, it was found that the shaded segments in figure 1 were governed by the mini- 
mum gage constraint. Subsequently, all nodes of the three shaded segments were 
assigned the same design variable and the rest of the design variables were distributed 
among the other three segments. 

An initial configuration having uniform 0.127-cm-thick (0.05 in.) cover panels was 
used to start the design process both in SWIFT and WIDOWAC. Stress and flutter results 
for this initial configuration are given in table I. Because of the difference in the analyti- 
cal models, the maximum stress under a uniform load in WIDOWAC is 613 MN/m 2 
(89.0 ksi) compared with 690 MN/m 2 (100 ksi) in SWIFT. A more refined finite element 
model having 162 degrees of freedom was run with NASTRAN (ref. 22) and gave a maxi- 
mum stress of 645 MN/m 2 (93.5 ksi). Flutter speed was 12 percent lower in WIDOWAC 
than in SWIFT. Because of these differences, the WIDOWAC model is more nearly 
critical in flutter and the SWIFT model is more nearly critical in stress. 

Table II shows the results of the optimization with seven design variables of the 
same initial design by the two programs. * Results are shown for the strength-only 
optimum design, the flutter-only optimum design, and the combined strength -flutter 
optimum design. Running times for WIDOWAC are shorter in spite of the extra com- 
plexity of the wing model. (A faster version of SWIFT is now available (as used in ref. 23) 
which takes only 295 CPU seconds for the combined strength-flutter design.) Because of 
the small number (seven) of design variables used in the present example the run times 
with Newton's method are only lower by about 50 percent than run times with the DFP 
method. 


In WIDOWAC, it was not possible to model the full-depth sandwich wing up to the 
leading and trailing edges where the thickness is zero. A narrow strip was therefore cut 
from both edges. The mass of the two strips, based on minimum gage thickness (com- 
patible with SWIFT results), is 104 kg (230 lb) and was added to WIDOWAC results. 
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Built-Up Wing 


A second program which was used for comparison with WIDOW AC is DAWNS 
(ref. 24). DAWNS (Design of Aircraft WiNg Structures) can be used for the fully stressed 
design of wings under stress, buckling, and minimum gage constraints. It generates its 
own loads for a symmetric pull-up maneuver using Mach box supersonic aerodynamics. 

A built-up delta wing with 65° sweep optimized by DAWNS (ref. 25) was used for compari- 
son with WIDOWAC. The wing planform and the final distribution of the skin thickness is 
shown in figure 3; the broken line surrounds the part of the wing and carry-through struc- 
ture that was modeled for the structural idealization. The finite element model used in 
DAWNS was modified by replacing each quadrilateral membrane element by two triangular 
elements. The applied loads were taken from the DAWNS run where they were calculated 
for a 3g maneuver at Mach 3. Many of the spars and ribs are critical in buckling; hence, 
buckling stress rather than yield stress governs their design. WIDOWAC does not calcu- 
late buckling stresses or loads; therefore, to obtain a meaningful comparison with DAWNS, 
the buckling stresses calculated by DAWNS were used as maximum stress constraints 
whenever they were lower than the yield stress. The wing was modeled by 86 triangular 
skin elements and 101 shear web elements totaling 156 degrees of freedom and was 
designed under stress and minimum gage constraints. A few runs with 6 to 10 design 
variables were made to find the critical areas of the wing and then the wing was optimized 
with 23 design variables allocated according to the information obtained by the earlier 
runs. Using the DFP algorithm, about 2300 analyses were required to converge the design 
which took 2400 seconds of CDC 6600 CPU time; rim time with Newton's method was only 
about 400 seconds for the 333 analyses required for convergence. Sixteen of the design 
variables were allocated to the skin panels and 7 to the shear webs. The final mass was 
9550 kg (21030 lbm) compared with a final mass of 8080 kg (17 800 lbm) in DAWNS. The 
skin thickness distribution was similar as is shown in figures 3 and 4. To investigate the 
reason for the big differences in mass, the DAWNS final design was analyzed by using 
WIDOWAC. Although the difference between the maximum displacement was only 2.5 per- 
cent, the maximum stress in WIDOWAC was 28 percent higher. These higher stresses 
made the WIDOWAC design heavier. The difference in the stresses is attributed to the 
different finite elements used in the two programs. DAWNS used a quadrilateral linear 
strain membrane element while WIDOWAC uses a triangular constant strain membrane 
element. The use of quadrilateral elements for this wing is not appropriate because some 
of these quadrilaterals are not planar. As shown in reference 26, the error induced by the 
use of nonplanar elements can be large. (The wing analyzed in ref. 26 is the same one as 
used in this paper.) 
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Number of Design Variables 


As was noted before, it is impractical to assign each finite element a design vari- 
able when using mathematical programing methods for the optimization. The technique 
used in this work for reducing the number of design variables is to divide the cover panels 
into segments and assign design variables to corners. A study was conducted to find the 
dependence of the final mass on the number of design variables. It should be noted, how- 
ever, that when the number of the design variables is smaller than that of the finite ele- 
ments there are many ways to assign the design variables. (As a result the curves shown 
in figs. 5 and 6 are not unique.) 

The first study was on the full-depth sandwich wing divided into segments as shown 
in figure 1. Figure 5 shows the dependence of the final mass on the number of design 
variables. The division into segments was held constant. Because the shaded segments 
in figure 1 always came out as minimum gage they were lumped together and assigned one 
design variable. The other three segments could have a maximum of 11 design variables. 
It is seen from figure 5 that for the given segmentation there is no need to use the full 
number of design variables. The result of designating the thickness of each finite element 
as a design variable was an additional 5 percent reduction in mass as shown by the mini- 
mum mass line in figure 5. This 5 percent represents the constraint of freezing the seg- 
mentation. It should be noted that, as the number of design variables was increased, the 
changes in the final design had the character of refinement rather than changes in design 
concept. 

The second study was on the built-up wing optimized for stress constraints where 
both the cover panels and the rib and spars are optimized. The strategy in assigning 
design variables to the rib and spar elements was according to the level of stress. Ele- 
ments with similar stress margins (u cr /o) were grouped together. As the number of 
design variables was increased the stress margins of previous designs were used. 

Figure 6 shows the dependence of the mass of the built-up wing on the number of 
design variables. Of the 23 design variables used for the last point on the curve, 16 were 
cover panel design variables and 7 were rib and spar design variables. The final design 
with 23 design variables as shown in figure 3 is not much different from the DAWNS final, 
fully stressed design with 160 design variables (one for each element). Some of the dif- 
ference may be attributed to the different finite elements used in DAWNS and WIDOW AC. 

From these two studies, it appears that it may be possible to get a good approxima- 
tion to the optimum design with a number of design variables which is much smaller than 
the number of finite elements. The two studies were done on relatively crude finite ele- 
ment models, and for more refined models, the ratio between the number of finite ele- 
ments and the number of design variables can be bigger. Because working with a rela- 
tively small number of design variables does involve a penalty in the final mass, these 


16 



procedures for reducing the number of design variables may be satisfactory for prelimi- 
nary design but not for the final design process. 


Discontinuity of Flutter Speed 

It is known (ref. 27) that the flutter speed may be a discontinuous function of air 
density or Mach number. A similar discontinuous behavior as a function of structural 
properties was revealed during an investigation of the utility of a high-modulus — low- 
density material (beryllium) being used in wing design. In this supplemental study, seg- 
ment III of the full-depth sandwich wing (fig. 1) was assumed to be made of beryllium and 
the rest of the wing was assumed to be titanium. Although no benefit was found from this 
use of beryllium (possibly because buckling was not considered), it was found that the 
flutter speed is discontinuous as a function of the patch thickness. 

Figure 7(a) shows the dependence of the lowest flutter speed on the beryllium patch 
thickness. The rest of the cover panels are uniform 0.254-cm-thick (0.1 in.) titanium. 
The discontinuity is the result of the appearance of a new flutter mode when the patch 
thickness exceeds 0.371 cm (0.15 in.). The flutter frequency and the first three natural 
frequencies are shown in figure 7(b). The flutter eigenvector for both modes is a combi- 
nation of primarily the first three natural vibration modes. For a patch thickness of 
0.42 cm (0.165 in.) and unit amplitude of a 6-mode flutter vector the amplitudes of the 
first three modes in the higher speed mode are 

aj = 0.594 a 2 = 0.431 exp(-0.411i) a 3 = 0.677 exp(-0.304i) 

The amplitudes of the lower speed mode are 

a^ = 0.673 a 2 = 0.280 exp(-2.22i) a 3 = 0.684 exp(-0.468i) 

The amplitudes of the fourth and higher natural vibration modes are small. Note that all 
three vibration modes participate significantly in both flutter modes. Actual V-g plots 
were not calculated, but on the basis of past experience, the switch between the two flutter 
modes may have taken place as shown in figure 8. 

In the present work, the flutter constraint is defined in terms of the flutter speed, 
and, therefore, the discontinuity in the flutter speed causes the breakdown of the optimi- 
zation procedure; the flutter constraint, however, may be formulated in terms of continu- 
ous parameters of the flutter phenomenon. For example, instead of equation (2) the fol- 
lowing expression may be used: 
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( 22 ) 


where g cr is the critical damping value and V m ^ n is a lower bound on the flutter 
speed. Other continuous constraints are possible (ref. 28); however, such constraints 
require more computation than the constraint of equation (2). 


CONCLUSIONS 


The present work explored methods for reducing the amount of computation involved 
in the automated design of complex structures for stress and flutter requirements. A 
pilot finite element program for the minimum mass design of wing structures under 
stress, flutter, and minimum gage constraints was developed. The following conclusions 
were reached: 

1. The use of appropriate analysis algorithms can reduce the computation effort 
required for a reanalysis of a structure to a small fraction of the computation effort 
required for the analysis of the original structure. 

2. When constraints are represented by a penalty function, Newton's method with 
approximate second derivatives may be much more efficient than the commonly used 
quasi-Newton methods (such as Davidon-Fletcher-Powell algorithm). 

3. The number of design variables can be significantly less than the number of 
finite elements without excessive mass penalties. 

4. The flutter speed may not be a continuous function of structural stiffness. There- 
fore, in optimization under flutter constraint, the constraint should be formulated in terms 
of other parameters of the flutter phenomenon which are continuous. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Hampton, Va., May 22, 1973. 
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APPENDIX A 


CONVERSION OF U.S. CUSTOMARY UNITS TO SI UNITS 

The International System of Units (SI) was adopted by the Eleventh General Con- 
ference on Weights and Measures, Paris, October 1960 (ref. 16). Conversion factors for 
the units herein are given in the following table: 


Physical quantity 

U.S. Customary 
Unit 

Conversion 

factor 

(*) 

SI Unit 

Force 

lbf 

4.44822 

newtons (N) 


f ft 

0.3048 

1 , v 

Length 

( in. 

0.0254 

>meters (m) 

Mass 

lbm 

0.45359 

kilograms (kg) 

Stress, modulus . . . 

ksi 

6.895 X 106 

newtons/meter 2 (N/m 2 ) 

Density 

lbm/in 2 

27.68 x 10 6 

grams/meter 2 (g/m 2 ) 

Pressure 

lbf/ft 2 

47.88 

newtons/meter 2 (N/m2) 

Speed 

ft/sec 

0.3048 

meters/second (m/s) 


^Multiply value given in U.S. Customary Unit by conversion factor to obtain 
equivalent value in SI Unit. 


Prefixes to indicate multiples of units are as follows: 


Prefix 

Multiple 

giga (G) 

109 

mega (M) 

106 

kilo (k) 

10 2 

deci (d) 

10' 1 

centi (c) 

lO- 2 

milli (m) 

lO- 2 
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APPENDIX B 


FINITE ELEMENT ANALYSIS OF SYMMETRIC WING 


WIDOW AC assumes that the wing has a symmetric profile and symmetric behavior; 
that is, the lower and upper surfaces of the wing are assumed to have x,y displacements 
of equal magnitude but opposite sign and the same z displacements. Therefore, only 
nodes on the upper surface have to be considered. The shear web element in the pro- 
gram is derived from the corresponding NASTRAN element (ref. 22, p. 5.3-1). It is 
defined by two nodes on the upper surface of the wing. The other two nodes are assumed 
to be on the z = 0 surface, to be tied rigidly to the upper surface nodes in the z direc- 
tion, and to be constrained from any movement in the x,y directions. 

If the upper two nodes have the coordinates (x^,y^,z^, (x2,y2, z 2)> then the stiffness 
matrix for the element [Kg] of thickness t e is given by 

[Ke] = ^ (eg) (cgf (Bl) 


where 


Z = 


( Z 1 + z 2) L 

Gt e 



2 

3(1 + v) 



(B2) 


and 


{/-g} ~ 2 A y L ’ ”( z l +z 2)’ x x L > A y^j z l +z 2 


IT 


L = (x 2 - x l} 2 + (y 2 - yi ) 2 ] 


1/2 


(B3) 


X x = 


x 2 - X 1 



y2 - y_i 

L 


J 


The shear flow is given by 

<ii * I ( c s) T ( w ) 


(B4) 
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where {w} is the displacement vector containing the three displacements at point 1 
followed by the three displacements at point 2. The shear stress s in the element is 
given as 

i = Y (B5) 
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TABLE I.- COMPARISON OF RESULTS OBTAINED FROM SWIFT AND WIDOW AC 
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CPU = CDC 6600 Central Processor time 










Constraints : 

1. Maximum stress 862 MN/m^ (125 ksi) 

2 

for uniform loading 6.89 kN/m (1 psi) 

2. Flutter speed at 7620 m (25 000 ft) 
not less than 762 m/s (2500 ft/sec) 

3. Minimum skin thickness, 0.0381 cm 

(0.015 in.) 

Figure 2.- Definition of full-depth sandwich wing and design constraints. 
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Thickness 
0.076 to 0.102 



Figure 3.- Cover panel thickness distribution obtained with WIDOW AC for built-up 65° delta wing 









Thickness 

0.076 to 0.102 cr 
(0.03 to 0.04 in.) 



Figure 4.- Cover-panel thickness distribution obtained with DAWNS for built-up 65° delta wing (ref. 25) 
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(a) Flutter speed. 
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Figure 7.- Natural frequency, flutter frequency, and flutter speed as a function 
of patch thickness for titanium -beryllium sandwich wing. 





Figure 8.- A sketch showing discontinuity of flutter 
speed on a V-g diagram. 


NASA-Langley, 1973 32 L-8592 



NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 
WASHINGTON, D.C. 20546 


OFFICIAL BUSINESS 

PENALTY FOR PRIVATE USE $300 SPECIAL FOURTH-CLASS RATE 

BOOK 


POSTAGE AND FEES PA 
NATIONAL AERONAUTICS 
SPACE ADMINISTRATE 
451 



POSTMASTER : 


If Undeliverable (Section 158 
Postal Manual) Do Not Return 


"The aeronautical and space activities of the United States shall be 
conducted so. as to contribute ... to the expansion of human knowl- 
edge of phenomena in the atmosphere and space. The Administration 
shall provide for the widest practicable and appropriate dissemination 
of information concerning its activities and the results thereof.” 

—National Aeronautics and Space Act of 1958 


NASA SCIENTIFIC AND TECHNICAL PUBLICATIONS 


TECHNICAL REPORTS: Scientific and 
technical information considered important, 
complete, and a lasting contribution to existing 
knowledge. 

TECHNICAL NOTES: Information less broad 
in scope but nevertheless of importance as a 
contribution to existing knowledge. 

TECHNICAL MEMORANDUMS: 
Information receiving limited distribution 
because of preliminary data, security classifica- 
tion, or other reasons. Also includes conference 
proceedings with either limited or unlimited 
distribution. 

CONTRACTOR REPORTS: Scientific and 
technical information generated under a NASA 
contract or grant and considered an important 
contribution to existing knowledge. 


TECHNICAL TRANSLATIONS: Information 
published in a foreign language considered 
to merit NASA distribution in English. 

SPECIAL PUBLICATIONS: Information 
derived from or of value to NASA activities. 
Publications include final reports of major 
projects, monographs, data compilations, 
handbooks, sourcebooks, and special 
bibliographies. 

TECHNOLOGY UTILIZATION 
PUBLICATIONS: Information on technology 
used by NASA that may be of particular 
interest in commercial and other non-aerospace 
applications. Publications include Tech Briefs, 
Technology Utilization Reports and 
Technology Surveys. 


Details on the availability of these publications may be obtained from: 

SCIENTIFIC AND TECHNICAL INFORMATION OFFICE 
NATIONAL AERONAUTICS AND SPACE ADMINISTRATION 

Washington, D.C. 20546 




